{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numba import cuda, jit, int32, float32, int64\n",
    "from numba.cuda.random import create_xoroshiro128p_states, xoroshiro128p_uniform_float32\n",
    "from math import pow, hypot, ceil\n",
    "import random\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Read The problem data file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "class vrp():\n",
    "    def __init__(self, capacity=None):\n",
    "        self.capacity = capacity\n",
    "        self.nodes = np.zeros((1,4), dtype=np.float32)\n",
    "    def addNode(self, label, demand, posX, posY):\n",
    "        newrow = np.array([label, demand, posX, posY], dtype=np.float32)\n",
    "        self.nodes = np.vstack((self.nodes, newrow))\n",
    "\n",
    "# Read the problem data file\n",
    "def readInput():\n",
    "\t# Create VRP object:\n",
    "    vrpManager = vrp()\n",
    "\t## First reading the VRP from the input ##\n",
    "    print('Reading data file...', end=' ')\n",
    "    fo = open('/home/conda_user/GA_VRP/test_set/P/P-n16-k8.vrp',\"r\")\n",
    "    lines = fo.readlines()\n",
    "    for i, line in enumerate(lines):\n",
    "        while line.upper().startswith('CAPACITY'):\n",
    "            inputs = line.split()\n",
    "            vrpManager.capacity = np.float32(inputs[2])\n",
    "\t\t\t# Validating positive non-zero capacity\n",
    "            if vrpManager.capacity <= 0:\n",
    "                print(sys.stderr, 'Invalid input: capacity must be neither negative nor zero!')\n",
    "                exit(1)\n",
    "            break       \n",
    "        while line.upper().startswith('NODE_COORD_SECTION'):\n",
    "            i += 1\n",
    "            line = lines[i]\n",
    "            while not (line.upper().startswith('DEMAND_SECTION') or line=='\\n'):\n",
    "                inputs = line.split()\n",
    "                vrpManager.addNode(np.int16(inputs[0]), 0.0, np.float32(inputs[1]), np.float32((inputs[2])))\n",
    "                # print(vrpManager.nodes)\n",
    "                i += 1\n",
    "                line = lines[i]\n",
    "                while (line=='\\n'):\n",
    "                    i += 1\n",
    "                    line = lines[i]\n",
    "                    if line.upper().startswith('DEMAND_SECTION'): break \n",
    "                if line.upper().startswith('DEMAND_SECTION'):\n",
    "                    i += 1\n",
    "                    line = lines[i] \n",
    "                    while not (line.upper().startswith('DEPOT_SECTION')):                  \n",
    "                        inputs = line.split()\n",
    "\t\t\t\t\t\t# Validating demand not greater than capacity\n",
    "                        if float(inputs[1]) > vrpManager.capacity:\n",
    "                            print(sys.stderr,\n",
    "\t\t\t\t\t\t\t'Invalid input: the demand of the node %s is greater than the vehicle capacity!' % vrpManager.nodes[0])\n",
    "                            exit(1)\n",
    "                        if float(inputs[1]) < 0:\n",
    "                            print(sys.stderr,\n",
    "                            'Invalid input: the demand of the node %s cannot be negative!' % vrpManager.nodes[0])\n",
    "                            exit(1)                            \n",
    "                        vrpManager.nodes[int(inputs[0])][1] =  float(inputs[1])\n",
    "                        i += 1\n",
    "                        line = lines[i]\n",
    "                        while (line=='\\n'):\n",
    "                            i += 1\n",
    "                            line = lines[i]\n",
    "                            if line.upper().startswith('DEPOT_SECTION'): break\n",
    "                        if line.upper().startswith('DEPOT_SECTION'):\n",
    "                            vrpManager.nodes = np.delete(vrpManager.nodes, 0, 0)                          \n",
    "                            print('Done.')\n",
    "                            return(vrpManager.capacity, vrpManager.nodes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Calculate fitness:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### CPU version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# define fitness kernel here:\n",
    "# @jit(nopython=True)\n",
    "def fitness(cost_table, individual):\n",
    "    zero_arr = np.zeros(1, dtype=np.int32)\n",
    "    zeroed_indiv = np.copy(individual)\n",
    "    \n",
    "    # nodes represent the row/column index in the cost table\n",
    "    for i in range(len(zeroed_indiv)):\n",
    "        zeroed_indiv[i] = zeroed_indiv[i] - 1\n",
    "        \n",
    "    if zeroed_indiv[0] != 0:\n",
    "        zeroed_indiv = np.hstack((zero_arr, zeroed_indiv))\n",
    "    if individual[-1] != 1:\n",
    "        zeroed_indiv = np.hstack((zeroed_indiv, zero_arr))\n",
    "        \n",
    "    fitness_val = 0\n",
    "    for i in range(len(zeroed_indiv)-1):\n",
    "        fitness_val += cost_table[int(zeroed_indiv[i]), int(zeroed_indiv[i+1])]\n",
    "        \n",
    "    return(fitness_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### GPU version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# define fitness kernel here:\n",
    "@cuda.jit\n",
    "def fitness_gpu(cost_table_d, individual_d, zeroed_indiv_d, fitness_val_d):     \n",
    "    \n",
    "    # nodes represent the row/column index in the cost table\n",
    "    threadId_row, threadId_col = cuda.grid(2)\n",
    "    \n",
    "    # Mapping between the 2D-grid indexing and the 1D-vector indexing:\n",
    "    index = threadId_row*(cuda.blockDim.x)+threadId_col\n",
    "    \n",
    "    fitness_val_d[0] = 0\n",
    "    if index+1 <= len(individual_d):\n",
    "        zeroed_indiv_d[index] = individual_d[index] - 1\n",
    "    \n",
    "    if index == 0 and zeroed_indiv_d[index] != 0:\n",
    "        cuda.atomic.add(fitness_val_d,0,cost_table_d[0, zeroed_indiv_d[index]])\n",
    "        cuda.atomic.add(fitness_val_d,0,cost_table_d[zeroed_indiv_d[index], zeroed_indiv_d[index+1]])\n",
    "    elif index == len(zeroed_indiv_d)-1 and zeroed_indiv_d[index] != 0:\n",
    "        cuda.atomic.add(fitness_val_d,0,cost_table_d[zeroed_indiv_d[index], 0])\n",
    "    elif index == len(zeroed_indiv_d)-1 and zeroed_indiv_d[index] == 0:\n",
    "        pass\n",
    "    elif index+1 <= len(zeroed_indiv_d):\n",
    "        cuda.atomic.add(fitness_val_d,0,cost_table_d[zeroed_indiv_d[index], zeroed_indiv_d[index+1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Every individual MUST be initialized with length 2 * no._of_nodes\n",
    "\n",
    "individual = np.array([8,15,1,5,12,1,14,9,1,11,16,16,1,6,6,4,1,2,1,3,1,7,1,\\\n",
    "                      8,15,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], dtype=np.int32)\n",
    "# individual = np.array(range(10000), dtype=np.int32)\n",
    "individual_d = cuda.to_device(individual)\n",
    "zeroed_indiv_d = cuda.to_device(individual)\n",
    "\n",
    "fitness_val_d = cuda.to_device(np.array([0], dtype=np.int32))\n",
    "\n",
    "fitness_gpu[blocks,threads_per_block](cost_table_d, individual_d, zeroed_indiv_d, fitness_val_d)\n",
    "print(fitness_val_d.copy_to_host()[0])\n",
    "\n",
    "###############################################################################################\n",
    "# Speed test of CPU and GPU versions of the function:\n",
    "# print(\"CPU time:\")\n",
    "# cost_table = np.zeros((data.shape[0], data.shape[0]), dtype=np.int32)\n",
    "# cost_table = calc_cost(data, popsize, vrp_capacity, cost_table)\n",
    "# %timeit fitness(cost_table, individual)\n",
    "# print(\"GPU time:\")\n",
    "# %timeit fitness_gpu[blocks,threads_per_block](cost_table_d, individual_d, zeroed_indiv_d, fitness_val_d)\n",
    "###############################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate cost table:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### CPU version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "## Calculate cost table:\n",
    "# @jit(nopython=True)\n",
    "def calc_cost(data, popsize, vrp_capacity, cost_table):\n",
    "    shifted_data = np.copy(data)\n",
    "    for i in range(len(shifted_data[:,0])):\n",
    "        shifted_data[i,0] = shifted_data[i,0] - 1\n",
    "\n",
    "    for row in range(len(shifted_data[:,0])):\n",
    "        for col in range(len(shifted_data[:,0])):\n",
    "            cost_table[row, col] = round(hypot((shifted_data[row,2] - shifted_data[col,2]),\\\n",
    "                                                       (shifted_data[row,3] - shifted_data[col,3])))\n",
    "    return cost_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GPU version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate cost table:\n",
    "@cuda.jit\n",
    "def calc_cost_gpu(data_d, popsize, vrp_capacity, cost_table_d):\n",
    "    threadId_row, threadId_col = cuda.grid(2)\n",
    "    \n",
    "#     data_d[threadId_row,0] = data_d[threadId_row,0] - 1\n",
    "    \n",
    "####ceil() is used instead of round() as the latter crashes the kernel.\n",
    "####This causes +1 values in some cost distances\n",
    "\n",
    "    if (threadId_row <= data_d.shape[0]-1) and (threadId_col <= data_d.shape[0]-1):\n",
    "        cost_table_d[threadId_row, threadId_col] = ceil(hypot(data_d[threadId_row,2] - data_d[threadId_col,2],\\\n",
    "                                                              data_d[threadId_row,3] - data_d[threadId_col,3]))\n",
    "#     popArr = initializePop(data, popsize, vrp_capacity, cost_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### CPU version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# @jit(nopython=True)\n",
    "def adjust(individual, data, vrp_capacity, cost_table):\n",
    "\n",
    "    # Delete duplicate nodes\n",
    "    adjusted_indiv = np.zeros(individual.shape[0], dtype=np.int32)\n",
    "    j = 0\n",
    "    for i in range(len(individual)):\n",
    "        if not np.any(individual[i] == adjusted_indiv):\n",
    "            adjusted_indiv[j] = individual[i]\n",
    "            j += 1\n",
    "\n",
    "    # Delete ones and zeros\n",
    "    adjusted_indiv = np.delete(adjusted_indiv, np.where(adjusted_indiv==1)[0])\n",
    "    adjusted_indiv = np.delete(adjusted_indiv, np.where(adjusted_indiv==0)[0])\n",
    "\n",
    "    \n",
    "    # Insert missing nodes\n",
    "    for i in range(data.shape[0]):\n",
    "        if not np.any(data[i,0] == adjusted_indiv):\n",
    "            adjusted_indiv = np.hstack((adjusted_indiv, np.array([data[i,0]], dtype=np.int32)))\n",
    "\n",
    "    i = 0               # index\n",
    "    reqcap = 0.0        # required capacity\n",
    "\n",
    "    while i < len(adjusted_indiv): \n",
    "        if adjusted_indiv[i] != 1:\n",
    "            reqcap += data[data[:,0] == adjusted_indiv[i]][0,1]\n",
    "        else:\n",
    "            reqcap = 0\n",
    "        \n",
    "        if reqcap > vrp_capacity: \n",
    "            adjusted_indiv = np.hstack((adjusted_indiv[:i], np.array([1], dtype=np.int32), adjusted_indiv[i:]))\n",
    "            reqcap = 0.0\n",
    "        i += 1\n",
    "        \n",
    "    if adjusted_indiv[0] != 1:\n",
    "        adjusted_indiv = np.hstack((np.array([1], dtype=np.int32), adjusted_indiv))\n",
    "    if adjusted_indiv[-1] != 1:\n",
    "        adjusted_indiv = np.hstack((adjusted_indiv, np.array([1], dtype=np.int32)))\n",
    "    \n",
    "#     adjusted_indiv = np.hstack((adjusted_indiv, np.asarray([fitness(cost_table, adjusted_indiv)], dtype=np.int32)))\n",
    "    return adjusted_indiv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### GPU version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "@cuda.jit\n",
    "def adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d, pop_d, zeroed_indiv_d, fitness_val_d):\n",
    "    \n",
    "    # nodes represent the row/column index in the cost table\n",
    "    threadId_row, threadId_col = cuda.grid(2)\n",
    "    \n",
    "    # Remove duplicated elements from every single individual/row in population array:\n",
    "    r_flag = 9999 # A flag for removal/replacement\n",
    "    \n",
    "    if threadId_row <= pop_d.shape[0]-1 and threadId_col <= pop_d.shape[1]-1 and threadId_col != 0:\n",
    "                    \n",
    "        for i in range(threadId_col-1, -1, -1):\n",
    "            if pop_d[threadId_row, threadId_col] == pop_d[threadId_row, i]\\\n",
    "            and pop_d[threadId_row, threadId_col] != 0:\n",
    "                pop_d[threadId_row, threadId_col] = r_flag \n",
    "            \n",
    "        for j in range(data_d.shape[0]):\n",
    "            for i in range(threadId_col-1, -1, -1):\n",
    "                if data_d[j,0] == pop_d[threadId_row, i]:\n",
    "                    missing_d[threadId_row, j] = 0\n",
    "                    break\n",
    "                else:\n",
    "                    missing_d[threadId_row, j] = data_d[j,0]\n",
    "                     \n",
    "    # Add missing nodes to every single individual:\n",
    "            \n",
    "    if threadId_col == pop_d.shape[1]-1:\n",
    "        missing_elements = True\n",
    "        for i in range(missing_d.shape[1]):\n",
    "                if missing_d[threadId_row, i] != 0:\n",
    "                    missing_elements = True\n",
    "                    for j in range(pop_d.shape[1]):\n",
    "                        if pop_d[threadId_row, j] == r_flag:\n",
    "                            pop_d[threadId_row, j] = missing_d[threadId_row, i]\n",
    "                            missing_d[threadId_row, i] = 0\n",
    "                            break\n",
    "                else:\n",
    "                    missing_elements = False\n",
    "\n",
    "        if not missing_elements:\n",
    "        # shift individual's elements to the left for every inserted '1':\n",
    "            for i in range(pop_d.shape[1], 0, -1):\n",
    "                if pop_d[threadId_row, i] == r_flag:\n",
    "                    for j in range(i, pop_d.shape[1]-1):\n",
    "                        new_val = pop_d[threadId_row, j+1]\n",
    "                        pop_d[threadId_row, j] = new_val\n",
    "\n",
    "        reqcap = 0.0        # required capacity\n",
    "        for i in range(pop_d.shape[1]-1):\n",
    "            if pop_d[threadId_row, i] != 1 and pop_d[threadId_row, i] != 0:\n",
    "                reqcap += data_d[pop_d[threadId_row, i]-1, 1]\n",
    "                if reqcap > vrp_capacity:\n",
    "                    # here will be the insert '1' algorithm:\n",
    "                    new_val = 1\n",
    "                    rep_val = pop_d[threadId_row, i]\n",
    "                    \n",
    "                    # shift individual's elements to the right for every inserted '1': \n",
    "                    for j in range(i, pop_d.shape[1]-1):\n",
    "                        pop_d[threadId_row, j] = new_val\n",
    "                        new_val = rep_val\n",
    "                        rep_val = pop_d[threadId_row, j+1]\n",
    "                    reqcap = 0.0                    \n",
    "            else:\n",
    "                reqcap = 0.0\n",
    "                \n",
    "            \n",
    "#         # The last part is to add the individual's fitness value at the very end of it.\n",
    "#             individual = pop_d[threadId_row,:]\n",
    "#             fitness_gpu(cost_table_d, individual, zeroed_indiv_d, fitness_val_d)            \n",
    "#             pop_d[threadId_row, -1] = fitness_val_d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### CPU version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def evolvePop(pop, vrp_data, iterations, popsize, vrp_capacity, extended_cost, opt, cost_table=0):\n",
    "    # Running the genetic algorithm\n",
    "    run_time = timer()\n",
    "    stucking_indicator = 0\n",
    "    for i in tqdm(range(iterations)):\n",
    "        old_best = pop[0][-1]\n",
    "        nextPop = []\n",
    "        nextPop_set = set()\n",
    "\n",
    "        elite_count = len(pop)//20      \n",
    "        sorted_pop = pop.copy()\n",
    "\n",
    "        # Apply two-opt for the new top 5% individuals:\n",
    "        for idx, individual in enumerate(sorted_pop[:elite_count]):\n",
    "            if individual[0] >= i:\n",
    "                sorted_pop[idx], cost = two_opt.two_opt(individual[1:-1], cost_table)\n",
    "                sorted_pop[idx].append(9999)\n",
    "                fitness_value = fitness(cost_table, sorted_pop[idx][:-1])\n",
    "                sorted_pop[idx][-1] = (fitness_value)\n",
    "                sorted_pop[idx].insert(0,individual[0])\n",
    "        \n",
    "        sorted_pop.sort(key= lambda elem: elem[-1])\n",
    "        pop = sorted_pop.copy()\n",
    "        \n",
    "        start_evolution_timer = timer()\n",
    "        # terminate if optimal is reached or runtime exceeds 1h\n",
    "        if ((sorted_pop[0][-1] + extended_cost) > opt) and (timer() - run_time <= 60):\n",
    "            nextPop = sorted_pop[:elite_count] # top 5% of the parents will remain in the new generation         \n",
    "\n",
    "            # for j in range(round(((len(pop))-elite_count) / 2)):\n",
    "            while len(nextPop_set) < popsize:\n",
    "                # Selecting randomly 4 individuals to select 2 parents by a binary tournament\n",
    "                parentIds = set()\n",
    "                while len(parentIds) < 4:\n",
    "                    parentIds.add(random.randint(0, len(pop) - 1))\n",
    "\n",
    "                # Avoid stucking to a local minimum swap after 25 generations of identical fitness\n",
    "                #if stucking_indicator >= 25:\n",
    "                    #print('\\nstucking is spotted', pop[1])\n",
    "                    #for idx, swapped_indiv in enumerate(pop[1:elite_count]):\n",
    "                        #i1 = swapped_indiv[1:round(len(swapped_indiv)/2)]\n",
    "                        #i2 = swapped_indiv[round(len(swapped_indiv)/2): -1]\n",
    "                        ## i1 = random.randint(1, len(swapped_indiv) - 2)\n",
    "                        ## i2 = random.randint(1, len(swapped_indiv) - 2)\n",
    "                        #swapped_indiv = i2\n",
    "                        #swapped_indiv = np.append(swapped_indiv, i1)\n",
    "                        ## swapped_indiv[i1], swapped_indiv[i2] = swapped_indiv[i2], swapped_indiv[i1]\n",
    "                        #swapped_indiv = adjust(np.asarray(swapped_indiv[1:], dtype=np.float32), np.asarray(vrp_data, dtype=np.float32), vrp_capacity)\n",
    "                        #fitness_val = fitness(np.asarray(vrp_data, np.float32), np.asarray(swapped_indiv[1:], np.float32))\n",
    "                        ## swapped_indiv[-1] = fitness_val\n",
    "                        #swapped_indiv = np.append(swapped_indiv, fitness_val)\n",
    "                        #pop[idx] = swapped_indiv\n",
    "                    #stucking_indicator = 0\n",
    "               \n",
    "                parentIds = list(parentIds)\n",
    "                # Selecting 2 parents with the binary tournament\n",
    "                parent1 = list(pop[parentIds[0]] if pop[parentIds[0]][len(pop[parentIds[0]])-1] < pop[parentIds[1]][len(pop[parentIds[1]])-1] else pop[parentIds[1]])\n",
    "                parent2 = list(pop[parentIds[2]] if pop[parentIds[2]][len(pop[parentIds[2]])-1] < pop[parentIds[3]][len(pop[parentIds[3]])-1] else pop[parentIds[3]])\n",
    "\n",
    "                child1 = parent1[1:].copy()\n",
    "                child2 = parent2[1:].copy()\n",
    "\n",
    "                # Performing Two-Point crossover and generating two children\n",
    "                # Selecting (n/5 - 1) random cutting points for crossover, with the same points (indexes) for both parents, based on the shortest parent\n",
    "\n",
    "                cutIdx = [0] * ((min(len(parent1) - 2, len(parent2) - 2))//5 - 1)\n",
    "                for k in range(0, len(cutIdx)):\n",
    "                    cutIdx[k] = random.randint(1, min(len(parent1) - 2, len(parent2) - 2))\n",
    "                    while cutIdx[k] in cutIdx[:k]:\n",
    "                        cutIdx[k] = random.randint(1, min(len(parent1) - 2, len(parent2) - 2))\n",
    "                cutIdx.sort()\n",
    "                for k in range(0, len(cutIdx), 2):\n",
    "                    if len(cutIdx) %2 == 1 and k == len(cutIdx) - 1: # Odd number\n",
    "                        child1[cutIdx[k]:] = child2[cutIdx[k]:]\n",
    "                        child2[cutIdx[k]:] = child1[cutIdx[k]:]\n",
    "                    else:                       \n",
    "                        child1[cutIdx[k]:cutIdx[k + 1]] = child2[cutIdx[k]:cutIdx[k + 1]]\n",
    "                        child2[cutIdx[k]:cutIdx[k + 1]] = child1[cutIdx[k]:cutIdx[k + 1]]        \n",
    "\n",
    "                # Doing mutation: swapping two positions in one of the individuals, with 1:15 probability\n",
    "                mutation_prob = 40\n",
    "                if random.randint(1, mutation_prob) == 1:\n",
    "                    # Random swap mutation\n",
    "                    ptomutate = child1\n",
    "                    i1 = random.randint(0, len(ptomutate) - 2)\n",
    "                    i2 = random.randint(0, len(ptomutate) - 2)\n",
    "                    # Repeat random selection if depot was selected\n",
    "                    while ptomutate[i1] == 1:\n",
    "                        i1 = random.randint(0, len(ptomutate) - 2)\n",
    "                    while ptomutate[i2] == 1:\n",
    "                        i2 = random.randint(0, len(ptomutate) - 2)\n",
    "                    ptomutate[i1], ptomutate[i2] = ptomutate[i2], ptomutate[i1]\n",
    "\n",
    "                if random.randint(1, mutation_prob) == 1:\n",
    "                    ptomutate = child2\n",
    "                    i1 = random.randint(0, len(ptomutate) - 2)\n",
    "                    i2 = random.randint(0, len(ptomutate) - 2)\n",
    "                    # Repeat random selection if depot was selected\n",
    "                    while ptomutate[i1] == 1:\n",
    "                        i1 = random.randint(0, len(ptomutate) - 2)\n",
    "                    while ptomutate[i2] == 1:\n",
    "                        i2 = random.randint(0, len(ptomutate) - 2)\n",
    "                    ptomutate[i1], ptomutate[i2] = ptomutate[i2], ptomutate[i1]\n",
    "\n",
    "                # Adjusting individuals               \n",
    "                child1 = adjust(np.asarray(child1, dtype=np.float32), np.asarray(vrp_data, dtype=np.float32), vrp_capacity)\n",
    "                child2 = adjust(np.asarray(child2, dtype=np.float32), np.asarray(vrp_data, dtype=np.float32), vrp_capacity)\n",
    "\n",
    "                # # Apply 2-opt:\n",
    "                # child1, fitness_val = two_opt.two_opt(child1[:-1], cost_table)\n",
    "                # child2, fitness_val = two_opt.two_opt(child2[:-1], cost_table)\n",
    "\n",
    "                fitness_val = fitness(cost_table, child1[:-1])\n",
    "                child1[-1] = fitness_val\n",
    "                \n",
    "                fitness_val = fitness(cost_table, child2[:-1])\n",
    "                child2[-1] = fitness_val\n",
    "\n",
    "                child1 = list(child1)\n",
    "                child2 = list(child2)\n",
    "\n",
    "                child1.insert(0, i + 1)\n",
    "                child2.insert(0, i + 1)\n",
    "\n",
    "                # Add children to population iff they are better than parents\n",
    "                if (child1[-1] < parent1[-1]) | (child1[-1] < parent2[-1]) | ((timer() - start_evolution_timer) > 30):\n",
    "                    nextPop_set.add(tuple(child1))\n",
    "                    # start_evolution_timer = timer()\n",
    "                    # nextPop_set.add(tuple(parent1))\n",
    "                \n",
    "                if (child2[-1] < parent1[-1]) | (child2[-1] < parent2[-1]) | ((timer() - start_evolution_timer) > 30):\n",
    "                    nextPop_set.add(tuple(child2))\n",
    "                    # start_evolution_timer = timer()\n",
    "                    # nextPop_set.add(tuple(parent2))   \n",
    "                               \n",
    "            nextPop = list(nextPop_set)\n",
    "\n",
    "            # Updating population generation\n",
    "\n",
    "            # random.shuffle(nextPop)\n",
    "            nextPop = sorted(nextPop, key= lambda elem: elem[-1])\n",
    "\n",
    "            if nextPop[0][-1] == old_best:\n",
    "                stucking_indicator += 1\n",
    "            else:\n",
    "                stucking_indicator = 0\n",
    "\n",
    "            pop = nextPop\n",
    "            if not (i+1) % 5: # print population every 300 generations\n",
    "                print(f'Population at generation {i+1}:{pop}\\nBest: {pop[0][-1]}')\n",
    "        elif (timer() - run_time >= 60):\n",
    "            print('Time criteria is met')\n",
    "            break\n",
    "        elif (((sorted_pop[0][-1] + extended_cost) <= opt)):\n",
    "            print('Cost criteria is met')\n",
    "            break\n",
    "    return (pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GPU version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit(device=True)\n",
    "def cut_points(pool_d, row, threadId_col):\n",
    "    \n",
    "    if threadId_col > 0 and threadId_col <= pool_d[row + 2, 2]:\n",
    "        pool_d[row + 3, threadId_col] = 10 # should be replaced by a random number\n",
    "        pool_d[row + 3, threadId_col+1] = 7\n",
    "        pool_d[row + 3, threadId_col+2] = 5\n",
    "\n",
    "# Crossover points (i.e., pool_d[row + 3,:]) must be ordered:\n",
    "    if threadId_col == 1:       \n",
    "        for i in range(1, pool_d[row + 2, 2]):\n",
    "            min_val = pool_d[row + 3, i]\n",
    "            min_index = i\n",
    "            for j in range(i, pool_d[row + 2, 2]+1):\n",
    "                if min_val > pool_d[row + 3, j]:\n",
    "                    min_val = pool_d[row + 3, j]\n",
    "                    min_index = j\n",
    "\n",
    "            pool_d[row + 3, i], pool_d[row + 3, min_index] = pool_d[row + 3, min_index], pool_d[row + 3, i]\n",
    "    cuda.syncthreads()\n",
    "# ----------------------------------------------------------------------------------------------------------\n",
    "@cuda.jit(device=True)\n",
    "def swap_indices(pool_d, row, threadId_col, j):\n",
    "    pool_d[row, threadId_col], pool_d[row+1, threadId_col] =\\\n",
    "    pool_d[row+1, threadId_col], pool_d[row, threadId_col]\n",
    "# ----------------------------------------------------------------------------------------------------------\n",
    "@cuda.jit(device=True)    \n",
    "def cross_over(generations, popsize, opt, vrp_capacity, cost_table_d, cut_idx_d, pool_d, new_pop_d, pop_d, row):\n",
    "    threadId_row, threadId_col = cuda.grid(2)\n",
    "    \n",
    "    if threadId_row < pop_d.shape[0]-1 and threadId_col <= pop_d.shape[1]-1 and threadId_col != 0:\n",
    "    #   Create a pool of 4 randomly selected individuals:\n",
    "        pool_d[row + 0, threadId_col] = pop_d[0, threadId_col] # 0 should be replaced by a random number\n",
    "        pool_d[row + 1, threadId_col] = pop_d[1, threadId_col] # 1 should be replaced by a random number\n",
    "        pool_d[row + 2, threadId_col] = pop_d[2, threadId_col] # 2 should be replaced by a random number\n",
    "        pool_d[row + 3, threadId_col] = pop_d[3, threadId_col] # 3 should be replaced by a random number\n",
    "\n",
    "\n",
    "    # Selecting 2 parents with the binary tournament\n",
    "    # The first two pool_d rows are re-assigned as parents:\n",
    "    # ----------------------------1st Parent--------------------------------------------------\n",
    "        if pool_d[row, -1] < pool_d[row + 1, -1]:\n",
    "            pool_d[row, threadId_col] = pool_d[row, threadId_col]\n",
    "        else:\n",
    "            pool_d[row, threadId_col] = pool_d[row + 1, threadId_col]\n",
    "\n",
    "    # ----------------------------2nd Parent--------------------------------------------------\n",
    "        if pool_d[row + 2, -1] < pool_d[row + 3, -1]:\n",
    "            pool_d[row + 1, threadId_col] = pool_d[row + 2, threadId_col]\n",
    "        else:\n",
    "            pool_d[row + 1, threadId_col] = pool_d[row + 3, threadId_col]\n",
    "\n",
    "        pool_d[row + 2, threadId_col] = 0\n",
    "        pool_d[row + 3, threadId_col] = 0\n",
    "\n",
    "    # Performing Two-Point crossover and generating two children:\n",
    "    # Calculate the actual length of parents, put it in pool_d\n",
    "        if pool_d[row, threadId_col] != 0:\n",
    "            cuda.atomic.add(pool_d, (row + 2, pool_d.shape[1]-2), 1)\n",
    "\n",
    "        if pool_d[row + 1, threadId_col] != 0:\n",
    "            cuda.atomic.add(pool_d, (row + 3, pool_d.shape[1]-2), 1)\n",
    "\n",
    "        pool_d[row + 2, 1] = \\\n",
    "        min(pool_d[row+2, pool_d.shape[1]-2], pool_d[row+3, pool_d.shape[1]-2]) # Minimum length of the two parents\n",
    "\n",
    "    # Select (n/5 - 1) random cutting points for crossover based on the shortest parent\n",
    "        pool_d[row + 2, 2] = pool_d[row + 2, 1]//5 - 1 # number of cutting points\n",
    "        cut_points(pool_d, row, threadId_col)\n",
    "# For odd swap indices, swap the chromosomes with indices less than the swap value:\n",
    "\n",
    "#     if pool_d[row + 2, 2] %2 == 1: # Number of indices is odd\n",
    "\n",
    "    if pool_d[row + 2, 2]%2 == 1:\n",
    "        for j in range(1, (pool_d[row + 2, 2])+1):\n",
    "            if threadId_col <= pool_d[row + 3, j] and threadId_col > pool_d[row + 3, j-1] and j%2 == 1:    \n",
    "                swap_indices(pool_d, row, threadId_col, j)\n",
    "    else:\n",
    "        for j in range(1, (pool_d[row + 2, 2])+1):\n",
    "            if (threadId_col <= pool_d[row + 3, j] and threadId_col > pool_d[row + 3, j-1] and j%2 == 1) or\\\n",
    "            (threadId_col > pool_d[row + 3, j] and j == pool_d[row + 2, 2]):    \n",
    "                swap_indices(pool_d, row, threadId_col, j)\n",
    "# ----------------------------------------------------------------------------------------------------------\n",
    "@cuda.jit(device=True)\n",
    "def mutate(pool_d, row, threadId_col):\n",
    "# Mutation: swapping two positions in the children, with 1:40 probability\n",
    "    mutation_prob = 40\n",
    "#     if random.randint(1, mutation_prob) == 1:\n",
    "    \n",
    "    if threadId_col == 1:    \n",
    "# Repeat random selection if depot was selected:    \n",
    "        i1 = 1\n",
    "        while pool_d[row, i1] == 1:\n",
    "            i1 = 4\n",
    "    #         i1 = random.randint(0, len(ptomutate) - 2)\n",
    "\n",
    "        i2 = 1\n",
    "        while pool_d[row, i2] == 1:\n",
    "            i2 = 15\n",
    "    #         i2 = random.randint(0, len(ptomutate) - 2)\n",
    "\n",
    "\n",
    "        pool_d[row, i1], pool_d[row, i2] = pool_d[row, i2], pool_d[row, i1]\n",
    "\n",
    "    # Repeat for the second child:    \n",
    "        i1 = 1\n",
    "        while pool_d[row+1, i1] == 1:\n",
    "            i1 = 4\n",
    "    #         i1 = random.randint(0, len(ptomutate) - 2)\n",
    "\n",
    "        i2 = 1\n",
    "        while pool_d[row+1, i2] == 1:\n",
    "            i2 = 16\n",
    "    #         i2 = random.randint(0, len(ptomutate) - 2)\n",
    "\n",
    "        pool_d[row+1, i1], pool_d[row+1, i2] = pool_d[row+1, i2], pool_d[row+1, i1]\n",
    "        cuda.syncthreads()\n",
    "# ----------------------------------------------------------------------------------------------------------\n",
    "\n",
    "@cuda.jit\n",
    "def evolvePop_gpu(count, generations, popsize, opt, vrp_capacity, data_d, cost_table_d,\\\n",
    "                  cut_idx_d, pool_d, new_pop_d, pop_d, fitness_val_d):\n",
    "    \n",
    "    # nodes represent the row/column index in the cost table\n",
    "    threadId_row, threadId_col = cuda.grid(2)\n",
    "    stride = 4\n",
    "    row = threadId_row*stride    \n",
    "    \n",
    "    if threadId_row < pop_d.shape[0]-1 and threadId_col <= pop_d.shape[1]-1 and threadId_col != 0:\n",
    "        cross_over(generations, popsize, opt, vrp_capacity, cost_table_d, cut_idx_d, pool_d,\\\n",
    "                   new_pop_d, pop_d, row)\n",
    "        mutate(pool_d, row, threadId_col)\n",
    "        adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d, pop_d, \\\n",
    "                   zeroed_indiv_d, fitness_val_d, threadId_row, threadId_col)\n",
    "        \n",
    "        new_pop_d[threadId_row, 0] = count\n",
    "        new_pop_d[threadId_row, threadId_col] = pool_d[row, threadId_col]\n",
    "\n",
    "#     # Running the genetic algorithm\n",
    "#     run_time = timer()\n",
    "#     stucking_indicator = 0\n",
    "#     for i in tqdm(range(iterations)):\n",
    "#         old_best = pop[0][-1]\n",
    "#         nextPop = []\n",
    "#         nextPop_set = set()\n",
    "\n",
    "#         elite_count = len(pop)//20      \n",
    "#         sorted_pop = pop.copy()\n",
    "\n",
    "# # Apply two-opt for the new top 5% individuals:\n",
    "#         for idx, individual in enumerate(sorted_pop[:elite_count]):\n",
    "#             if individual[0] >= i:\n",
    "#                 sorted_pop[idx], cost = two_opt.two_opt(individual[1:-1], cost_table)\n",
    "#                 sorted_pop[idx].append(9999)\n",
    "#                 fitness_value = fitness(cost_table, sorted_pop[idx][:-1])\n",
    "#                 sorted_pop[idx][-1] = (fitness_value)\n",
    "#                 sorted_pop[idx].insert(0,individual[0])\n",
    "        \n",
    "#         sorted_pop.sort(key= lambda elem: elem[-1])\n",
    "#         pop = sorted_pop.copy()\n",
    "        \n",
    "#         start_evolution_timer = timer()\n",
    "# # terminate if optimal is reached or runtime exceeds 1h\n",
    "#         if ((sorted_pop[0][-1] + extended_cost) > opt) and (timer() - run_time <= 60):\n",
    "#             nextPop = sorted_pop[:elite_count] # top 5% of the parents will remain in the new generation         \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# # Adjusting individuals               \n",
    "#                 child1 = adjust(np.asarray(child1, dtype=np.float32), np.asarray(vrp_data, dtype=np.float32), vrp_capacity)\n",
    "#                 child2 = adjust(np.asarray(child2, dtype=np.float32), np.asarray(vrp_data, dtype=np.float32), vrp_capacity)\n",
    "\n",
    "#                 fitness_val = fitness(cost_table, child1[:-1])\n",
    "#                 child1[-1] = fitness_val\n",
    "                \n",
    "#                 fitness_val = fitness(cost_table, child2[:-1])\n",
    "#                 child2[-1] = fitness_val\n",
    "\n",
    "#                 child1 = list(child1)\n",
    "#                 child2 = list(child2)\n",
    "\n",
    "# # Insert generation number at the beginning of every individual:\n",
    "#                 child1.insert(0, i + 1)\n",
    "#                 child2.insert(0, i + 1)\n",
    "\n",
    "#                 # Add children to population iff they are better than parents\n",
    "#                 if (child1[-1] < parent1[-1]) | (child1[-1] < parent2[-1]) | ((timer() - start_evolution_timer) > 30):\n",
    "#                     nextPop_set.add(tuple(child1))\n",
    "#                     # start_evolution_timer = timer()\n",
    "#                     # nextPop_set.add(tuple(parent1))\n",
    "                \n",
    "#                 if (child2[-1] < parent1[-1]) | (child2[-1] < parent2[-1]) | ((timer() - start_evolution_timer) > 30):\n",
    "#                     nextPop_set.add(tuple(child2))\n",
    "#                     # start_evolution_timer = timer()\n",
    "#                     # nextPop_set.add(tuple(parent2))   \n",
    "                               \n",
    "#             nextPop = list(nextPop_set)\n",
    "\n",
    "#             # Updating population generation\n",
    "\n",
    "#             # random.shuffle(nextPop)\n",
    "#             nextPop = sorted(nextPop, key= lambda elem: elem[-1])\n",
    "\n",
    "#             if nextPop[0][-1] == old_best:\n",
    "#                 stucking_indicator += 1\n",
    "#             else:\n",
    "#                 stucking_indicator = 0\n",
    "\n",
    "#             pop = nextPop\n",
    "#             if not (i+1) % 5: # print population every 300 generations\n",
    "#                 print(f'Population at generation {i+1}:{pop}\\nBest: {pop[0][-1]}')\n",
    "#         elif (timer() - run_time >= 60):\n",
    "#             print('Time criteria is met')\n",
    "#             break\n",
    "#         elif (((sorted_pop[0][-1] + extended_cost) <= opt)):\n",
    "#             print('Cost criteria is met')\n",
    "#             break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypingError",
     "evalue": "Failed in nopython mode pipeline (step: nopython frontend)\nInvalid use of Function(<numba.cuda.compiler.DeviceFunctionTemplate object at 0x7f0a530b1cc0>) with argument(s) of type(s): (array(float32, 2d, A), float32, array(int32, 2d, A), readonly array(int32, 2d, A), array(int32, 2d, A), readonly array(int32, 1d, A), array(int32, 2d, A), int32, int32)\n * parameterized\nIn definition 0:\n    TypingError: Failed in nopython mode pipeline (step: nopython frontend)\nInvalid use of Function(<built-in function setitem>) with argument(s) of type(s): (readonly array(int32, 2d, A), (int32, int64), Literal[int](0))\n * parameterized\nIn definition 0:\n    All templates rejected with literals.\nIn definition 1:\n    All templates rejected without literals.\nIn definition 2:\n    All templates rejected with literals.\nIn definition 3:\n    All templates rejected without literals.\nIn definition 4:\n    All templates rejected with literals.\nIn definition 5:\n    All templates rejected without literals.\nIn definition 6:\n    TypeError: Cannot modify value of type readonly array(int32, 2d, A)\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typing/arraydecl.py:179\nIn definition 7:\n    TypeError: Cannot modify value of type readonly array(int32, 2d, A)\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typing/arraydecl.py:179\nThis error is usually caused by passing an argument of a type that is unsupported by the named function.\n[1] During: typing of setitem at <ipython-input-117-8b0ca80f3242> (44)\n\nFile \"<ipython-input-117-8b0ca80f3242>\", line 44:\ndef adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d,\\\n    <source elided>\n                if data_d[j,0] == pop_d[threadId_row, i]:\n                    missing_d[threadId_row, j] = 0\n                    ^\n\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typeinfer.py:915\nIn definition 1:\n    TypingError: Failed in nopython mode pipeline (step: nopython frontend)\nInvalid use of Function(<built-in function setitem>) with argument(s) of type(s): (readonly array(int32, 2d, A), (int32, int64), Literal[int](0))\n * parameterized\nIn definition 0:\n    All templates rejected with literals.\nIn definition 1:\n    All templates rejected without literals.\nIn definition 2:\n    All templates rejected with literals.\nIn definition 3:\n    All templates rejected without literals.\nIn definition 4:\n    All templates rejected with literals.\nIn definition 5:\n    All templates rejected without literals.\nIn definition 6:\n    TypeError: Cannot modify value of type readonly array(int32, 2d, A)\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typing/arraydecl.py:179\nIn definition 7:\n    TypeError: Cannot modify value of type readonly array(int32, 2d, A)\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typing/arraydecl.py:179\nThis error is usually caused by passing an argument of a type that is unsupported by the named function.\n[1] During: typing of setitem at <ipython-input-117-8b0ca80f3242> (44)\n\nFile \"<ipython-input-117-8b0ca80f3242>\", line 44:\ndef adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d,\\\n    <source elided>\n                if data_d[j,0] == pop_d[threadId_row, i]:\n                    missing_d[threadId_row, j] = 0\n                    ^\n\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typeinfer.py:915\nThis error is usually caused by passing an argument of a type that is unsupported by the named function.\n[1] During: resolving callee type: Function(<numba.cuda.compiler.DeviceFunctionTemplate object at 0x7f0a530b1cc0>)\n[2] During: typing of call at <ipython-input-121-614642961486> (134)\n\n\nFile \"<ipython-input-121-614642961486>\", line 134:\ndef evolvePop_gpu(count, generations, popsize, opt, vrp_capacity, data_d, cost_table_d,\\\n    <source elided>\n        adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d, pop_d, \\\n                   zeroed_indiv_d, fitness_val_d, threadId_row, threadId_col)\n                   ^\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypingError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-122-1d789771adbd>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0mevolvePop_gpu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mblocks\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mthreads_per_block\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m              (count, generations, popsize, opt, vrp_capacity, data_d, cost_table_d, \\\n\u001b[0;32m---> 12\u001b[0;31m               cut_idx_d, pool_d, new_pop_d, pop_d, fitness_val_d)\n\u001b[0m\u001b[1;32m     13\u001b[0m     \u001b[0mcuda\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msynchronize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m     \u001b[0mcount\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/cuda/compiler.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args)\u001b[0m\n\u001b[1;32m    804\u001b[0m         \u001b[0mSpecialize\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0minvoke\u001b[0m \u001b[0mthis\u001b[0m \u001b[0mkernel\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    805\u001b[0m         '''\n\u001b[0;32m--> 806\u001b[0;31m         \u001b[0mkernel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mspecialize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    807\u001b[0m         \u001b[0mcfg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkernel\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgriddim\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mblockdim\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstream\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msharedmem\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    808\u001b[0m         \u001b[0mcfg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/cuda/compiler.py\u001b[0m in \u001b[0;36mspecialize\u001b[0;34m(self, *args)\u001b[0m\n\u001b[1;32m    815\u001b[0m         argtypes = tuple(\n\u001b[1;32m    816\u001b[0m             [self.typingctx.resolve_argument_type(a) for a in args])\n\u001b[0;32m--> 817\u001b[0;31m         \u001b[0mkernel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margtypes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    818\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mkernel\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    819\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/cuda/compiler.py\u001b[0m in \u001b[0;36mcompile\u001b[0;34m(self, sig)\u001b[0m\n\u001b[1;32m    831\u001b[0m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtargetoptions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'link'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    832\u001b[0m             kernel = compile_kernel(self.py_func, argtypes,\n\u001b[0;32m--> 833\u001b[0;31m                                     **self.targetoptions)\n\u001b[0m\u001b[1;32m    834\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdefinitions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margtypes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkernel\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    835\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbind\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler_lock.py\u001b[0m in \u001b[0;36m_acquire_compile_lock\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m     30\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0m_acquire_compile_lock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_acquire_compile_lock\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/cuda/compiler.py\u001b[0m in \u001b[0;36mcompile_kernel\u001b[0;34m(pyfunc, args, link, debug, inline, fastmath, extensions, max_registers)\u001b[0m\n\u001b[1;32m     60\u001b[0m def compile_kernel(pyfunc, args, link, debug=False, inline=False,\n\u001b[1;32m     61\u001b[0m                    fastmath=False, extensions=[], max_registers=None):\n\u001b[0;32m---> 62\u001b[0;31m     \u001b[0mcres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcompile_cuda\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpyfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtypes\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvoid\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdebug\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdebug\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minline\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minline\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     63\u001b[0m     \u001b[0mfname\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcres\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfndesc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mllvm_func_name\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     64\u001b[0m     lib, kernel = cres.target_context.prepare_cuda_kernel(cres.library, fname,\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler_lock.py\u001b[0m in \u001b[0;36m_acquire_compile_lock\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m     30\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0m_acquire_compile_lock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_acquire_compile_lock\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/cuda/compiler.py\u001b[0m in \u001b[0;36mcompile_cuda\u001b[0;34m(pyfunc, return_type, args, debug, inline)\u001b[0m\n\u001b[1;32m     49\u001b[0m                                   \u001b[0mreturn_type\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mreturn_type\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m                                   \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mflags\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 51\u001b[0;31m                                   locals={})\n\u001b[0m\u001b[1;32m     52\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     53\u001b[0m     \u001b[0mlibrary\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcres\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlibrary\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler.py\u001b[0m in \u001b[0;36mcompile_extra\u001b[0;34m(typingctx, targetctx, func, args, return_type, flags, locals, library, pipeline_class)\u001b[0m\n\u001b[1;32m    953\u001b[0m     pipeline = pipeline_class(typingctx, targetctx, library,\n\u001b[1;32m    954\u001b[0m                               args, return_type, flags, locals)\n\u001b[0;32m--> 955\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mpipeline\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompile_extra\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    956\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    957\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler.py\u001b[0m in \u001b[0;36mcompile_extra\u001b[0;34m(self, func)\u001b[0m\n\u001b[1;32m    375\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlifted\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    376\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlifted_from\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 377\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_compile_bytecode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    378\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    379\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mcompile_ir\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunc_ir\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlifted\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlifted_from\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler.py\u001b[0m in \u001b[0;36m_compile_bytecode\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    884\u001b[0m         \"\"\"\n\u001b[1;32m    885\u001b[0m         \u001b[0;32massert\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunc_ir\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 886\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_compile_core\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    887\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    888\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_compile_ir\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler.py\u001b[0m in \u001b[0;36m_compile_core\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    871\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdefine_pipelines\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    872\u001b[0m         \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfinalize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 873\u001b[0;31m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstatus\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    874\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mres\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    875\u001b[0m             \u001b[0;31m# Early pipeline completion\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler_lock.py\u001b[0m in \u001b[0;36m_acquire_compile_lock\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m     30\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0m_acquire_compile_lock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_acquire_compile_lock\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, status)\u001b[0m\n\u001b[1;32m    252\u001b[0m                     \u001b[0;31m# No more fallback pipelines?\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    253\u001b[0m                     \u001b[0;32mif\u001b[0m \u001b[0mis_final_pipeline\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 254\u001b[0;31m                         \u001b[0;32mraise\u001b[0m \u001b[0mpatched_exception\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    255\u001b[0m                     \u001b[0;31m# Go to next fallback pipeline\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    256\u001b[0m                     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, status)\u001b[0m\n\u001b[1;32m    243\u001b[0m                 \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    244\u001b[0m                     \u001b[0mevent\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"-- %s\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mstage_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 245\u001b[0;31m                     \u001b[0mstage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    246\u001b[0m                 \u001b[0;32mexcept\u001b[0m \u001b[0m_EarlyPipelineCompletion\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    247\u001b[0m                     \u001b[0;32mreturn\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler.py\u001b[0m in \u001b[0;36mstage_nopython_frontend\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    499\u001b[0m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    500\u001b[0m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturn_type\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 501\u001b[0;31m                 self.locals)\n\u001b[0m\u001b[1;32m    502\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtypemap\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtypemap\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    503\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturn_type\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreturn_type\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/compiler.py\u001b[0m in \u001b[0;36mtype_inference_stage\u001b[0;34m(typingctx, interp, args, return_type, locals)\u001b[0m\n\u001b[1;32m   1103\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1104\u001b[0m         \u001b[0minfer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_constraint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1105\u001b[0;31m         \u001b[0minfer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpropagate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1106\u001b[0m         \u001b[0mtypemap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrestype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcalltypes\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minfer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munify\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1107\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typeinfer.py\u001b[0m in \u001b[0;36mpropagate\u001b[0;34m(self, raise_errors)\u001b[0m\n\u001b[1;32m    913\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    914\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mraise_errors\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 915\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    916\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    917\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypingError\u001b[0m: Failed in nopython mode pipeline (step: nopython frontend)\nInvalid use of Function(<numba.cuda.compiler.DeviceFunctionTemplate object at 0x7f0a530b1cc0>) with argument(s) of type(s): (array(float32, 2d, A), float32, array(int32, 2d, A), readonly array(int32, 2d, A), array(int32, 2d, A), readonly array(int32, 1d, A), array(int32, 2d, A), int32, int32)\n * parameterized\nIn definition 0:\n    TypingError: Failed in nopython mode pipeline (step: nopython frontend)\nInvalid use of Function(<built-in function setitem>) with argument(s) of type(s): (readonly array(int32, 2d, A), (int32, int64), Literal[int](0))\n * parameterized\nIn definition 0:\n    All templates rejected with literals.\nIn definition 1:\n    All templates rejected without literals.\nIn definition 2:\n    All templates rejected with literals.\nIn definition 3:\n    All templates rejected without literals.\nIn definition 4:\n    All templates rejected with literals.\nIn definition 5:\n    All templates rejected without literals.\nIn definition 6:\n    TypeError: Cannot modify value of type readonly array(int32, 2d, A)\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typing/arraydecl.py:179\nIn definition 7:\n    TypeError: Cannot modify value of type readonly array(int32, 2d, A)\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typing/arraydecl.py:179\nThis error is usually caused by passing an argument of a type that is unsupported by the named function.\n[1] During: typing of setitem at <ipython-input-117-8b0ca80f3242> (44)\n\nFile \"<ipython-input-117-8b0ca80f3242>\", line 44:\ndef adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d,\\\n    <source elided>\n                if data_d[j,0] == pop_d[threadId_row, i]:\n                    missing_d[threadId_row, j] = 0\n                    ^\n\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typeinfer.py:915\nIn definition 1:\n    TypingError: Failed in nopython mode pipeline (step: nopython frontend)\nInvalid use of Function(<built-in function setitem>) with argument(s) of type(s): (readonly array(int32, 2d, A), (int32, int64), Literal[int](0))\n * parameterized\nIn definition 0:\n    All templates rejected with literals.\nIn definition 1:\n    All templates rejected without literals.\nIn definition 2:\n    All templates rejected with literals.\nIn definition 3:\n    All templates rejected without literals.\nIn definition 4:\n    All templates rejected with literals.\nIn definition 5:\n    All templates rejected without literals.\nIn definition 6:\n    TypeError: Cannot modify value of type readonly array(int32, 2d, A)\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typing/arraydecl.py:179\nIn definition 7:\n    TypeError: Cannot modify value of type readonly array(int32, 2d, A)\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typing/arraydecl.py:179\nThis error is usually caused by passing an argument of a type that is unsupported by the named function.\n[1] During: typing of setitem at <ipython-input-117-8b0ca80f3242> (44)\n\nFile \"<ipython-input-117-8b0ca80f3242>\", line 44:\ndef adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d,\\\n    <source elided>\n                if data_d[j,0] == pop_d[threadId_row, i]:\n                    missing_d[threadId_row, j] = 0\n                    ^\n\n    raised from /home/conda_user/anaconda3/envs/pycuda_env/lib/python3.7/site-packages/numba/typeinfer.py:915\nThis error is usually caused by passing an argument of a type that is unsupported by the named function.\n[1] During: resolving callee type: Function(<numba.cuda.compiler.DeviceFunctionTemplate object at 0x7f0a530b1cc0>)\n[2] During: typing of call at <ipython-input-121-614642961486> (134)\n\n\nFile \"<ipython-input-121-614642961486>\", line 134:\ndef evolvePop_gpu(count, generations, popsize, opt, vrp_capacity, data_d, cost_table_d,\\\n    <source elided>\n        adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d, pop_d, \\\n                   zeroed_indiv_d, fitness_val_d, threadId_row, threadId_col)\n                   ^\n"
     ]
    }
   ],
   "source": [
    "new_pop_d = cuda.device_array_like(pop_d)\n",
    "pool_d = np.zeros(shape=(5*popsize, pop_d.shape[1]), dtype=np.int32)\n",
    "pool_d = cuda.to_device(pool_d)\n",
    "cut_idx = [9999]*np.ones(shape=(pop_d.shape[1]), dtype=np.int32)\n",
    "cut_idx_d = cuda.to_device(cut_idx)\n",
    "\n",
    "count = 0\n",
    "\n",
    "while count <= generations:\n",
    "    evolvePop_gpu[blocks, threads_per_block]\\\n",
    "             (count, generations, popsize, opt, vrp_capacity, data_d, cost_table_d, \\\n",
    "              cut_idx_d, pool_d, new_pop_d, pop_d, fitness_val_d)\n",
    "    cuda.synchronize()\n",
    "    count += 1\n",
    "\n",
    "import sys\n",
    "np.set_printoptions(threshold=sys.maxsize)\n",
    "print(pool_d.copy_to_host()[52:56,:])\n",
    "print(new_pop_d.copy_to_host()[1:50,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adjust individuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit(device=True)\n",
    "def add_missing_nodes(r_flag, threadId_row, threadId_col, missing_d, pop_d):\n",
    "    # Add missing nodes to every single individual:\n",
    "            \n",
    "    if threadId_col == 1:\n",
    "        missing_elements = True\n",
    "        for i in range(missing_d.shape[1]):\n",
    "            if missing_d[threadId_row, i] != 0:\n",
    "                missing_elements = True\n",
    "                for j in range(pop_d.shape[1]):\n",
    "                    if pop_d[threadId_row, j] == r_flag:\n",
    "                        pop_d[threadId_row, j] = missing_d[threadId_row, i]\n",
    "                        missing_d[threadId_row, i] = 0\n",
    "                        break\n",
    "            else:\n",
    "                missing_elements = False\n",
    "\n",
    "#         if not missing_elements:\n",
    "#         # shift individual's elements to the left for every inserted '1':\n",
    "#             for i in range(pop_d.shape[1], 0, -1):\n",
    "#                 if pop_d[threadId_row, i] == r_flag:\n",
    "#                     for j in range(i, pop_d.shape[1]-1):\n",
    "#                         new_val = pop_d[threadId_row, j+1]\n",
    "#                         pop_d[threadId_row, j] = new_val\n",
    "        cuda.syncthreads()\n",
    "        \n",
    "@cuda.jit(device=True)\n",
    "def adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d,\\\n",
    "               pop_d, zeroed_indiv_d, fitness_val_d, threadId_row, threadId_col):\n",
    "    \n",
    "    # Remove duplicated elements from every single individual/row in population array:\n",
    "    r_flag = 9999 # A flag for removal/replacement\n",
    "    \n",
    "    if threadId_row <= pop_d.shape[0]-1 and threadId_col <= pop_d.shape[1]-1 and threadId_col != 0:\n",
    "                    \n",
    "        for i in range(threadId_col-1, -1, -1):\n",
    "            if pop_d[threadId_row, threadId_col] == pop_d[threadId_row, i]\\\n",
    "            and pop_d[threadId_row, threadId_col] != 0:\n",
    "                pop_d[threadId_row, threadId_col] = r_flag \n",
    "            \n",
    "        for j in range(data_d.shape[0]):\n",
    "            for i in range(threadId_col-1, -1, -1):\n",
    "                if data_d[j,0] == pop_d[threadId_row, i]:\n",
    "                    missing_d[threadId_row, j] = 0\n",
    "                    break\n",
    "                else:\n",
    "                    missing_d[threadId_row, j] = data_d[j,0]\n",
    "        \n",
    "        add_missing_nodes(r_flag, threadId_row, threadId_col, missing_d, pop_d)                     \n",
    "        cuda.syncthreads()\n",
    "\n",
    "#         reqcap = 0.0        # required capacity\n",
    "#         for i in range(pop_d.shape[1]-1):\n",
    "#             if pop_d[threadId_row, i] != 1 and pop_d[threadId_row, i] != 0:\n",
    "#                 reqcap += data_d[pop_d[threadId_row, i]-1, 1]\n",
    "#                 if reqcap > vrp_capacity:\n",
    "#                     # here will be the insert '1' algorithm:\n",
    "#                     new_val = 1\n",
    "#                     rep_val = pop_d[threadId_row, i]\n",
    "                    \n",
    "#                     # shift individual's elements to the right for every inserted '1': \n",
    "#                     for j in range(i, pop_d.shape[1]-1):\n",
    "#                         pop_d[threadId_row, j] = new_val\n",
    "#                         new_val = rep_val\n",
    "#                         rep_val = pop_d[threadId_row, j+1]\n",
    "#                     reqcap = 0.0                    \n",
    "#             else:\n",
    "#                 reqcap = 0.0\n",
    "                \n",
    "# --------------------------------------------------------------------------------------------            \n",
    "#         # The last part is to add the individual's fitness value at the very end of it.\n",
    "#             individual = pop_d[threadId_row,:]\n",
    "#             fitness_gpu(cost_table_d, individual, zeroed_indiv_d, fitness_val_d)            \n",
    "#             pop_d[threadId_row, -1] = fitness_val_d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialize population:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit\n",
    "def initializePop_gpu(rng_states, data_d, vrp_capacity, cost_table_d, missing_d,\\\n",
    "                      pop_d, zeroed_indiv_d, fitness_val_d):\n",
    "    \n",
    "    threadId_row, threadId_col = cuda.grid(2)\n",
    "    \n",
    "    # Generate the individuals from the nodes in data_d:\n",
    "    if threadId_col <= data_d.shape[0]-1:\n",
    "        pop_d[threadId_row, threadId_col] = data_d[threadId_col, 0]\n",
    "    \n",
    "    # Randmly shuffle each individual on a separate thread:   \n",
    "    col = 0\n",
    "    if threadId_row <= pop_d.shape[0]-1 and threadId_col <= data_d.shape[0]-1 and threadId_col != 0:\n",
    "        while col == 0:\n",
    "            rnd = (xoroshiro128p_uniform_float32(rng_states, threadId_row*threadId_col)*(data_d.shape[0]-1))\n",
    "            col = int(rnd)+1\n",
    "\n",
    "        pop_d[threadId_row, threadId_col], pop_d[threadId_row, col] =\\\n",
    "        pop_d[threadId_row, col], pop_d[threadId_row, threadId_col]\n",
    "        \n",
    "        adjust_gpu(data_d, vrp_capacity, cost_table_d, missing_d, pop_d,\\\n",
    "                   zeroed_indiv_d, fitness_val_d, threadId_row, threadId_col)\n",
    "        \n",
    "    # Adjust individuals using adjust_gpu function:\n",
    "    # Calculate fitness of each individual using fitness_gpu function:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading data file... Done.\n",
      "[[ 1  3 13  1  2  2  3 16  4  5  6  7  5  8  9 10  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]\n",
      " [ 1 15  1 14 11 10  5  4  3  6  2  3  4 16 12  9  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]\n",
      " [ 1 12  7  9 13  8  3 16  5  2  1 15  2  6  3  4  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]\n",
      " [ 1 16  5  6  3  4 12  1  2 14  8 15 13  2  3  4  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]\n",
      " [ 1 14  4 15 13 16 12  1  9  2  2 10  3 11  3  4  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]]\n",
      "-----------------------\n",
      "[[ 0  0  0  4  0  6  7  8  0  0 11 12  0 14 15  0  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  2  0  0  0  0  7  8  9  0  0 12 13  0  0  0  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  4  0  6  0  0  0 10 11  0  0 14  0  0  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  2  0  0  0  0  7  0  9 10 11  0  0  0  0  0  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  3  0  5  6  7  8  0  0 11  0  0  0  0  0  0  0  0  0  0  0  0  0\n",
      "   0  0  0  0  0  0  0  0  0  0]]\n"
     ]
    }
   ],
   "source": [
    "vrp_capacity, data = readInput()\n",
    "popsize = 100\n",
    "generations = 7000\n",
    "opt = 450\n",
    "\n",
    "data_d = cuda.to_device(data)\n",
    "cost_table_d = cuda.device_array(shape=(data.shape[0], data.shape[0]), dtype=np.int32)\n",
    "\n",
    "pop = np.zeros((popsize, 2*data.shape[0]+2), dtype=np.int32)\n",
    "pop_d = cuda.to_device(pop)\n",
    "\n",
    "zeros = np.zeros(shape=(popsize, pop_d.shape[1]), dtype=np.int32)\n",
    "missing_d = cuda.to_device(zeros)\n",
    "\n",
    "zeroed_indiv = np.zeros(shape=(pop_d.shape[1]), dtype=np.int32)\n",
    "zeroed_indiv_d = cuda.to_device(zeroed_indiv)\n",
    "fitness_val = np.zeros(shape=(popsize,1), dtype=np.int32)\n",
    "fitness_val_d = cuda.to_device(fitness_val)\n",
    "\n",
    "# GPU grid configurations:\n",
    "threads_per_block = (20, 20)\n",
    "blocks_no = ceil(max(popsize, 2*data.shape[0]+2)/threads_per_block[0])\n",
    "\n",
    "blocks = (blocks_no, blocks_no)\n",
    "rng_states = create_xoroshiro128p_states(threads_per_block[0]**2  * blocks[0]**2, seed=1)\n",
    "calc_cost_gpu[blocks, threads_per_block](data_d, popsize, vrp_capacity, cost_table_d)\n",
    "\n",
    "initializePop_gpu[blocks, threads_per_block]\\\n",
    "                 (rng_states, data_d, vrp_capacity, cost_table_d, missing_d,\\\n",
    "                  pop_d, zeroed_indiv_d, fitness_val_d)\n",
    "\n",
    "np.set_printoptions(threshold=sys.maxsize)\n",
    "print(pop_d.copy_to_host()[0:5,:], end='\\n-----------------------\\n')\n",
    "print(missing_d.copy_to_host()[0:5,:], end='\\n')\n",
    "\n",
    "# fitness_gpu[blocks,threads_per_block](cost_table_d, adjusted_indiv, zeroed_indiv_d, fitness_val_d)\n",
    "# print(fitness_val_d.copy_to_host()[0])\n",
    "# print(cost_table_d.copy_to_host())\n",
    "###############################################################################################\n",
    "# Speed test of CPU and GPU versions of the function:\n",
    "# cost_table = np.zeros((data.shape[0],data.shape[0]), dtype=np.int32)\n",
    "# print(calc_cost(data, popsize, vrp_capacity, cost_table).shape)\n",
    "# print('CPU time:')\n",
    "# %timeit calc_cost(data, popsize, vrp_capacity, cost_table)\n",
    "# print('GPU time:')\n",
    "#%timeit calc_cost_gpu[blocks, threads_per_block](data_d, popsize, vrp_capacity, cost_table_d)\n",
    "################################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# #%timeit adjust_gpu[blocks,threads_per_block]\\\n",
    "#  #(data_d, vrp_capacity, cost_table_d, missing_d, pop_d)\n",
    "# adjust_gpu[blocks, threads_per_block]\\\n",
    "# (data_d, vrp_capacity, cost_table_d, pop_d, zeroed_indiv_d, fitness_val_d)\n",
    "\n",
    "# # print(missing_d.copy_to_host()[50:65,:])\n",
    "# print(pop_d.copy_to_host()[0:4,:])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
